#RFS replication data and code###

#load dataset
#setwd("~/YOURFOLDER")

#load  libraries
library(cjpowR)
library(foreign)
library(readxl)
library(openxlsx)
library(stargazer)
library(tidyverse)
library(ggpubr)
library(psych)
library(car)
library(sjPlot)
library(panelr)
library(sjlabelled)
library(sjmisc)
library(plm)
library(gplots)
library(corrplot)
library(ggplot2)
library(ISLR)
library(ggrepel)
library(jsonlite)
library(pxR)
library(jtools)
library(viridis)
library(hrbrthemes)
library(BBmisc)
library(survey)
library(cjoint)
library(janitor)
library(reshape2)
library(data.table)
library(nnet)
library(xlsx)
library(readxl)
library(MASS)
library(ggpubr)
library(dplyr)
library(ltm)

#load data####
master <- readRDS(file = "RFSdata.rds")

#descriptives####
#table A1
age <- round(prop.table(table(master$group, master$Age),1),2)
gender <- round(prop.table(table(master$group, master$Gender),1),2)
edu <- round(prop.table(table(master$group, master$Education),1),2)
place <- round(prop.table(table(master$group, master$Place),1),2)
life <- round(prop.table(table(master$group, master$Life),1),2)
relig <- round(prop.table(table(master$group, master$Religion),1),2)
left <- round(prop.table(table(master$group, master$LeftRight2),1),2)
summaries <- as.data.frame(cbind(age,gender,edu,place,life,relig,left))
summaries$group <- as.factor(rownames(summaries))
write.xlsx(summaries, "group-summaries.xlsx")

#table A2
aff1 <- round(prop.table(table(master$group, master$Affected1_rec),1),2)
aff2 <- round(prop.table(table(master$group, master$Affected2_rec),1),2)
aff3 <- round(prop.table(table(master$group, master$Affected3_rec),1),2)
aff4 <- round(prop.table(table(master$group, master$Affected4_rec),1),2)
aff5 <- round(prop.table(table(master$group, master$Affected5_rec),1),2)
aff6 <- round(prop.table(table(master$group, master$Affected6_rec),1),2)
aff7 <- round(prop.table(table(master$group, master$Affected7_rec),1),2)
aff8 <- round(prop.table(table(master$group, master$Affected8_rec),1),2)
summaries <- as.data.frame(cbind(aff1,aff2,aff3,aff4,aff5,aff6,aff7,aff8))
summaries$group <- as.factor(rownames(summaries))
write.xlsx(summaries, "group-summariesaffected.xlsx")




#conjoint analyses####
Canada <- master[which(master$Country=='Canada'),]
Belgium <- master[which(master$Country=='Belgium'),]
France <- master[which(master$Country=='France'),]
Switzerland <- master[which(master$Country=='Switzerland' & (master$Language=='(Swiss-)German' |
                                                               master$Language=='French')),]
Australia <- master[which(master$Country=='Australia'),]
USA <- master[which(master$Country=='USA'),]


#bilingual countries
Canada$mothertongue <- NA
Canada$mothertongue[Canada$group=="CAN-French"] <- "French"
Canada$mothertongue[Canada$group=="CAN-English"] <- "English"
Canada$mothertongue <- factor(Canada$mothertongue,labels = c("French", "English"))
table(Canada$mothertongue) 

Belgium$mothertongue <- NA
Belgium$mothertongue[Belgium$group=="B-French"] <- "French"
Belgium$mothertongue[Belgium$group=="B-Dutch"] <- "Dutch"
Belgium$mothertongue <- factor(Belgium$mothertongue,labels = c("French", "Dutch"))
table(Belgium$mothertongue) 

France$mothertongue <- NA
France$mothertongue[France$Language=="French"] <- 1
France$mothertongue <- factor(France$mothertongue,labels = c("French"))
table(France$mothertongue) 

Switzerland$mothertongue <- NA
Switzerland$mothertongue[Switzerland$group=="CH-French"] <- "French"
Switzerland$mothertongue[Switzerland$group=="CH-German"] <- "German"
Switzerland$mothertongue <- as.factor(Switzerland$mothertongue)
table(Switzerland$mothertongue) 


#mono-language groups
francophones <- master[which(master$Language=='French' & !Country=="USA" & !Country=="Australia"),]
francophones$Country2 <- NA
francophones$Country2[francophones$Country=="Belgium"] <- "Belgium"
francophones$Country2[francophones$Country=="Switzerland"] <- "Switzerland"
francophones$Country2[francophones$Country=="Canada"] <- "Canada"
francophones$Country2[francophones$Country=="France"] <- "France"
francophones$Country2 <- as.factor(francophones$Country2)
table(francophones$Country2) 

#english-speakers
anglophones <- master[which(master$Language=='English'),]
anglophones$Country2 <- NA
anglophones$Country2[anglophones$Country=="Canada"] <- "Canada"
anglophones$Country2[anglophones$Country=="USA"] <- "USA"
anglophones$Country2[anglophones$Country=="Australia"] <- "Australia"
anglophones$Country2 <- as.factor(anglophones$Country2)
table(anglophones$Country2) 

#Figure 1
f1b <- cj(francophones, Y ~ govt + citizens + experts + stakeholders + transparency, id = ~ id, estimate = "mm", by= ~ Country2)
plot(f1b, group = "Country2", vline=0.5)+ theme_bw()+ggtitle("")+
  aes(shape = Country2) +
  scale_shape_manual(values=c(1, 2, 3, 4),na.translate = FALSE) +
  scale_colour_manual(values=c("black","red","green","blue"),na.translate = FALSE)+
  theme(legend.title = element_blank())


#Figure 2
f2 <- cj(anglophones, Y ~ govt + citizens + experts + stakeholders + transparency, id = ~ id, estimate = "mm", by= ~ Country2)
plot(f2, group = "Country2", vline=0.5)+ theme_bw()+ggtitle("")+
  aes(shape = Country2) +
  scale_shape_manual(values=c(1, 2, 3),na.translate = FALSE) +
  scale_colour_manual(values=c("black","red","green"),na.translate = FALSE)+
  theme(legend.title = element_blank())

#Figure 3
francophones$Country2 <- factor(francophones$Country2, levels = c("France", "Belgium", "Switzerland", "Canada"))

diff_mms <- cj(francophones, Y ~ govt + citizens + experts + stakeholders + transparency, id = ~ id, estimate = "mm_diff",
               by = ~Country2)

plot(diff_mms)+theme_bw()+ggtitle("")+ggplot2::facet_wrap(~BY, ncol = 3L)+
  scale_colour_manual(values=rep("black",5),na.translate = FALSE)+guides(color = "none")


#Figure 4
Canada$group5 <- NA
Canada$group5[Canada$Language=="English" & Canada$Country=="Canada"] <- "anglo-RoC"
Canada$group5[Canada$Language=="French" & Canada$Country=="Canada"] <- "franco-RoC"
Canada$group5[Canada$Language=="French" & Canada$Region=="Quebec"] <- "franco-quebecois"
Canada$group5[Canada$Language=="English" & Canada$Region=="Quebec"] <- "anglo-quebecois"
Canada$group5 <- as.factor(Canada$group5)

diff_mms <- cj(Canada, Y ~ govt + citizens + experts + stakeholders + transparency, id = ~ id, estimate = "mm_diff",
               by = ~group5)

plot(diff_mms)+theme_bw()+ggtitle("")+ggplot2::facet_wrap(~BY, ncol = 2L)+
  scale_colour_manual(values=rep("black",10),na.translate = FALSE)+guides(color = "none")

#Figure 5
#B
diff_mms <- cj(Belgium, Y ~ govt + citizens + experts + stakeholders + transparency, id = ~ id, estimate = "mm_diff",
               by = ~mothertongue)
plot(diff_mms)+theme_bw()+ggtitle("Belgium, 2 groups")+ggplot2::facet_wrap(~BY, ncol = 3L)+
  scale_colour_manual(values=rep("black",5),na.translate = FALSE)+guides(color = "none")+xlim(-.15,.15)

cj_anova(Belgium, Y ~ govt + citizens + experts + stakeholders + transparency, by= ~ mothertongue)

#CAN
diff_mms <- cj(Canada, Y ~ govt + citizens + experts + stakeholders + transparency, id = ~ id, estimate = "mm_diff",
               by = ~mothertongue)

plot(diff_mms)+theme_bw()+ggtitle("Canada, 2 groups")+ggplot2::facet_wrap(~BY, ncol = 3L)+
  scale_colour_manual(values=rep("black",5),na.translate = FALSE)+guides(color = "none")+xlim(-.15,.15)

cj_anova(Canada, Y ~ govt + citizens + experts + stakeholders + transparency, by= ~ mothertongue)


#CH
diff_mms <- cj(Switzerland, Y ~ govt + citizens + experts + stakeholders + transparency, id = ~ id, estimate = "mm_diff",
               by = ~mothertongue)
plot(diff_mms)+theme_bw()+ggtitle("Switzerland, 2 groups")+ggplot2::facet_wrap(~BY, ncol = 3L)+
  scale_colour_manual(values=rep("black",5),na.translate = FALSE)+guides(color = "none")

cj_anova(Switzerland, Y ~ govt + citizens + experts + stakeholders + transparency, by= ~ mothertongue)

#Table 3: F-statistics####
cj_anova(francophones, Y ~ govt + citizens + experts + stakeholders + transparency, by = ~ Country2)
cj_anova(francophones, Y ~ govt, by = ~ Country2)
cj_anova(francophones, Y ~ citizens, by = ~ Country2)
cj_anova(francophones, Y ~ experts, by = ~ Country2)
cj_anova(francophones, Y ~ stakeholders, by = ~ Country2)
cj_anova(francophones, Y ~ transparency, by = ~ Country2)

cj_anova(anglophones, Y ~ govt + citizens + experts + stakeholders + transparency, by = ~ Country2)
cj_anova(anglophones, Y ~ govt, by = ~ Country2)
cj_anova(anglophones, Y ~ citizens, by = ~ Country2)
cj_anova(anglophones, Y ~ experts, by = ~ Country2)
cj_anova(anglophones, Y ~ stakeholders, by = ~ Country2)
cj_anova(anglophones, Y ~ transparency, by = ~ Country2)

cj_anova(Belgium, Y ~ govt + citizens + experts + stakeholders + transparency, by = ~ mothertongue)
cj_anova(Belgium, Y ~ govt, by = ~ mothertongue)
cj_anova(Belgium, Y ~ citizens, by = ~ mothertongue)
cj_anova(Belgium, Y ~ experts, by = ~ mothertongue)
cj_anova(Belgium, Y ~ stakeholders, by = ~ mothertongue)
cj_anova(Belgium, Y ~ transparency, by = ~ mothertongue)

cj_anova(Canada, Y ~ govt + citizens + experts + stakeholders + transparency, by = ~ mothertongue)
cj_anova(Canada, Y ~ govt, by = ~ mothertongue)
cj_anova(Canada, Y ~ citizens, by = ~ mothertongue)
cj_anova(Canada, Y ~ experts, by = ~ mothertongue)
cj_anova(Canada, Y ~ stakeholders, by = ~ mothertongue)
cj_anova(Canada, Y ~ transparency, by = ~ mothertongue)

table(Switzerland$Language)
Switzerlandtwogroups <- subset(Switzerland, Language=="French" | Language=="(Swiss-)German")
cj_anova(Switzerlandtwogroups, Y ~ govt + citizens + experts + stakeholders + transparency, by = ~ mothertongue)
cj_anova(Switzerlandtwogroups, Y ~ govt, by = ~ mothertongue)
cj_anova(Switzerlandtwogroups, Y ~ citizens, by = ~ mothertongue)
cj_anova(Switzerlandtwogroups, Y ~ experts, by = ~ mothertongue)
cj_anova(Switzerlandtwogroups, Y ~ stakeholders, by = ~ mothertongue)
cj_anova(Switzerlandtwogroups, Y ~ transparency, by = ~ mothertongue)


#Figure A2####
f3 <- cj(master, Y ~ govt + citizens + experts + stakeholders + transparency, id = ~ id, estimate = "mm", by= ~ group)
plot(f3, group = "group", vline=0.5)+ theme_bw()+ggtitle("")+
  aes(shape = group) +
  scale_shape_manual(values=c(1, 2, 3,4,5,6,7,8,9,10),na.translate = FALSE) +
  scale_colour_manual(values=c("black",
                               "black","red",
                               "black","green",
                               "black","green",
                               "red","red"),na.translate = FALSE)+
  theme(legend.title = element_blank())

#Figure A3####
f3 <- cj(master, Y ~ govt + citizens + experts + stakeholders + transparency, id = ~ id, estimate = "mm", by= ~ group)
plot(f3, group = "group", vline=0.5)+ theme_bw()+ggtitle("")+
  aes(shape = group) +
  scale_shape_manual(values=c(1, 2, 3,4,5,6,7,8,9,10),na.translate = FALSE) +
  scale_colour_manual(values=c("black",
                               "black","red",
                               "black","green",
                               "black","green",
                               "red","red"),na.translate = FALSE)+
  theme(legend.position = "none")+facet_wrap(~BY, ncol = 5L)

#Figure A4, robustness 1: only matching respondents####
master$matching3 <- NA
master$matching3[master$matching=="French-Paris"] <- "F France"
master$matching3[master$matching=="French not Paris"] <- "F France"
master$matching3[master$matching=="Walloons"] <- "B Walloons"
master$matching3[master$matching=="Flemish"] <- "B Flemish"
master$matching3[master$matching=="Romands"] <- "CH Romands"
master$matching3[master$matching=="Deutschschweizer"] <- "CH Deutschschweizer"
master$matching3[master$matching=="Quebecois"] <- "CAN Quebecois"
master$matching3[master$matching=="Anglo-Canadians"] <- "CAN Anglo-Canadians"
master$matching3[master$Country=="Australia" & master$Language=="English"] <- "Australians"
master$matching3[master$Country=="USA" & master$Language=="English"] <- "US-Americans"

master$matching3 <- as.factor(master$matching3)

f4 <- cj(master, Y ~ govt + citizens + experts + stakeholders + transparency, id = ~ id, estimate = "mm", by= ~ matching3)
plot(f4, group = "matching3", vline=0.5)+ theme_bw()+ggtitle("")+
  aes(shape = matching3) +
  scale_shape_manual(values=c(1, 2, 3,4,5,6,7,8,9,10),na.translate = FALSE) +
  scale_colour_manual(values=c("red",
                               "green","black",
                               "red","black",
                               "green","black",
                               "black",
                               "red"),na.translate = FALSE)+
  theme(legend.position = "none")+facet_wrap(~BY, ncol = 5L)


#Figure A5: robustness 2, only passed covid test####
subset1 <- subset(master,CovidTest=="pass")

f5 <- cj(subset1, Y ~ govt + citizens + experts + stakeholders + transparency, id = ~ id, estimate = "mm", by= ~ group)
plot(f5, group = "group", vline=0.5)+ theme_bw()+ggtitle("")+
  aes(shape = group) +
  scale_shape_manual(values=c(1, 2, 3,4,5,6,7,8,9,10),na.translate = FALSE) +
  scale_colour_manual(values=c("black",
                               "black","red",
                               "black","green",
                               "black","green",
                               "red","red"),na.translate = FALSE)+
  theme(legend.position = "none")+facet_wrap(~BY, ncol = 5L)



#power analyses####
#Table A3
cjpowr_amcie(delta3 = 0.05, levels1=5, levels2=2, n=25984)
cjpowr_amcie(delta3 = 0.05, levels1=5, levels2=2, n=22712)
cjpowr_amcie(delta3 = 0.05, levels1=5, levels2=2, n=9760)
cjpowr_amcie(delta3 = 0.05, levels1=5, levels2=2, n=12096)
cjpowr_amcie(delta3 = 0.05, levels1=5, levels2=2, n=12160)

#Figure A5
d2 <- expand.grid(
  amce = seq(from=0.02, to = 0.1, by=0.001), 
  power = c(0.3, 0.5, 0.7, 0.9), 
  alpha = 0.05, 
  levels = 5,
  treat.prob = 0.5,
  sims = 10000)


df <- list2DF(do.call(cjpowr_amce, d2))
df$Power <- as.factor(df$power)
ggplot(df, aes(x = amce, y = n,color=Power,linetype=Power))+
  geom_line() +theme_bw()+
  xlab("AMCE") +
  ylab("Effective sample size")

